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ABSTRACT 

We demonstrate how cosmic star-formation history can be measured with one-point 
statistics of carbon-monoxide intensity maps. Using a P(D) analysis, the luminosity 
function of CO-emitting sources can be inferred from the measured one-point inten¬ 
sity PDF. The star-formation rate density (SFRD) can then be obtained, at several 
redshifts, from the CO luminosity density. We study the effects of instrumental noise, 
line foregrounds, and target redshift, and obtain constraints on the CO luminosity 
density of order 10%. We show that the SFRD uncertainty is dominated by that of 
the model connecting CO luminosity and star formation. For pessimistic estimates 
of this model uncertainty, we obtain an error of order 50% on SFRD for surveys 
targeting redshifts between 2 and 7 with reasonable noise and foregrounds included. 
However, comparisons between intensity maps and galaxies could substantially reduce 
this model uncertainty. In this case our constraints on SFRD at these redshifts improve 
to roughly 5 — 10%, which is highly competitive with current measurements. 

Key -words: cosmology: theory - galaxies: high-redshift 


1 INTRODUCTION 

In recent years, our understanding of the high-redshift uni¬ 
verse has dramatically improved thanks to a variety of multi¬ 
wavelength galaxy surveys. In particular, we have a much 
improved understanding of the history of star formation 
even in distant galaxies, indicating that the cosmic star for¬ 
mation rate density (SFRD) increases up to a maximum 
around 2 ~ 2 — 3 before declining until the present day 
(Madau & Dickinson (2014) and references therein). How¬ 
ever, most of our knowledge of very high redshifts comes 
from a relatively small population of the brightest galaxies. 

In this Letter we will discuss intensity mapping, a new 
means to study the distribution and properties of high- 
redshift galaxies. We will show how it is possible to obtain 
strong constraints on SFRD using this intensity mapping 
technique, in principle as strong as ~ 10% for near future 
surveys targeting 2 = 2.6. This constraint does not diminish 
significantly at higher redshifts, which suggests that inten¬ 
sity mapping may be a powerful tool for constraining the 
high-redshift SFRD. 

Instead of observing the emission of single galaxies, in¬ 
tensity mapping involves measuring how the intensity of a 
single spectral line varies on relatively large scales. With a 
beam much larger than a typical galaxy, an intensity map¬ 
ping experiment measures the aggregate emission from a 
large number of galaxies, requiring neither the high sensi¬ 
tivity or high resolution of a traditional galaxy survey. 

Intensity mapping is commonly discussed using the 21 
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cm neutral hydrogen line, specihcally targeted at the epoch 
of reionization (Madau et al. 1997). However in recent years 
observers have contemplated using other lines. Proposed 
lines include CH (Gong et al. 2012; Silva et al. 2014) and 
Lya (Pullen et al. 2014; Gong et al. 2014), but we will con¬ 
sider here the 115 GHz CO (1-0) rotational transition. CO 
is a commonly-observed molecule in nearby galaxies as it 
provides an excellent tracer of molecular gas, which in turn 
traces star formation activity (Bolatto et al. 2013). In an 
intensity mapping experiment, CO could be used to mea¬ 
sure star formation in high-redshift galaxies, without being 
restricted to the brightest sources. At least one such exper¬ 
iment, the Carbon Monoxide Mapping Array Pathfinder, is 
currently in the planning stages (Li et al. 2015). 

The typical statistic used to study an intensity map is 
the power spectrum P{k), which can be written as 

Pco{k, z) = (Tco)^ («)Pgal(fc, z) -f Pahot{z), (1) 

where (Tco) is the sky-averaged brightness temperature, 
Pgai is the galaxy power spectrum, and Pshot is scale- 
independent shot noise due to the discreteness of the CO 
sources. Both (Too) and Pshot depend sensitively on the 
properties of the galaxy population, but both depend only 
on integrals over the CO luminosity function, which limits 
their ability to constrain its shape. On top of that, (Too) is 
degenerate with the galaxy bias and the matter growth func¬ 
tion, and Pshot may be degenerate with shot noise contribu¬ 
tions from foreground lines such as HCN(l-O) (Breysse et al. 
2015). 

As the CO luminosity and star formation rate (SFR) 
of a galaxy may be related nonlinearly, these integral con¬ 
straints are suboptimal for constraining SFRD. What we 
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would like is a full measurement of the CO luminosity func¬ 
tion, which would allow us to make much better SFR esti¬ 
mates. For this, we turn to the one-point statistics of the 
map. Specifically, we use a technique known as P(D) anal¬ 
ysis, first proposed by Scheuer (1957). P(D) analysis is a 
method for obtaining source number counts in a confusion- 
limited survey. Some recent examples of its use include stud¬ 
ies of the cosmic infrared background (Glenn et al. 2010) 
and gamma-ray emission from our Galaxy (Lee et al. 2009, 
2015). Below we will demonstrate how this technique can 
relate the one-point PDF of an intensity map to the un¬ 
derlying luminosity function and illustrate its potential for 
constraining high-redshift star formation. 

2 P(D) ANALYSIS 

The one-point PDF of a CO intensity map depends both on 
the CO luminosity function and the clustering properties of 
the galaxy population. These quantities are related by 

OO 

P{T) = ^ P{N)Pm{T), (2) 

N=0 

where P{T) is the probability of observing a pixel with 
brightness temperature T, Pn{T) is the probability of ob¬ 
serving intensity T in a pixel containing N sources, and 
P{N) is the probability that a given pixel contains N 
sources (Lee et al. 2009). The PDF for a pixel containing 
zero sources is a delta function, and the PDF for a pixel 
with a single source can be determined from the luminosity 
function. For higher values of N, Pn{T) can be determined 
through a series of convolutions: 

POO 

Pn{T)= / Pn-i{T-T i)Pi{Ti)dTi. (3) 

Jo 

We then need a model for the luminosity function. We 
follow Breysse et al. (2015) and assume that the CO lumi¬ 
nosity L of a halo is related to its mass M by a power law, 
L = AM^, where Aco = 2x 10“® and hco = 1, and luminos¬ 
ity and mass are in Solar units. We consider only halos more 
massive than Mmin = 10® Solar masses, and assume that 
only a fraction /duty = tage/(10® yr) of them are emitting 
CO at any given time, where tage is the age of the universe in 
years. The number density of halos with a given mass is then 
determined by a mass function dn/dM (Tinker et al. 2008). 
If n is the mean number density of halos with masses greater 
than Mmin calculated by integrating this mass function, then 
we can write Pi(M) = {l/n)dn/dM. We can then convert 
this to Pi{L) using our mass-luminosity relation, and from 
there to a brightness temperature (Lidz et al. 2011). 

If we assume that the sources are unclustered, then 
the probability P{N) is a Poisson distribution with mean 
p = nPpix for pixels with volume Ppix. However, the source 
clustering will alter the shape of this distribution (Barcons 
1992). A perfect understanding of the clustering would re¬ 
quire knowledge of the full multi-point statistics of the 
source population. However we can obtain a reasonable es¬ 
timate of the clustering by assuming the sources are dis¬ 
tributed log-normally (Coles & Jones 1991). For this work, 
we follow the assumptions made in Breysse et al. (2015) and 
assume that the number of sources A in a pixel is deter¬ 
mined by drawing from a Poisson distribution Ppoiss(A, p) 
with mean p which is in turn drawn from a log-normal dis¬ 
tribution PhN{p)- We then have 


where 

PLN(/r) = / -[lnMn) + .l/2]P2.l ^ ^ 5 ^ 

/i^27r cr| 

and ffQ is the variance parameter from Equation (3.6) of 
Breysse et al. (2015), which depends on the chosen form of 
the power spectrum (however it is of order unity and has a 
small impact on the final form of P(T)). 

It should be noted that the proof-of-concept model con¬ 
sidered here is simplified in several respects for ease of com¬ 
putation. For example, the duty cycle /duty we use here is 
independent of halo mass, and we have a hard cutoff on 
SFR below Mmin. A more physical treatment would be to 
use a duty cycle that varies smoothly with mass and goes 
to zero at low masses (see for example Jaacks et al. (2012)). 
In addition, a realistic CO luminosity function likely has 
a “knee” at a fainter luminosity than our model predicts 
(Obreschkow et al. 2009), meaning that our calculation may 
over-predict the number counts of the brightest sources. It 
may be preferable to use a luminosity function model with 
more than two parameters which can accurately capture this 
“knee”, such as a Schechter function. We leave exploration 
of these issues and their effects on our results to future work 
(Breysse et ah, in preparation). 

With Equations (1-3) we can calculate the one-point 
PDF of a CO intensity map. However, a real map will include 
contributions from instrumental noise as well as foreground 
emission. For this Letter we ignore the effects of contin¬ 
uum foregrounds such as dust and synchrotron emission, but 
we do consider foregrounds with line spectra (Breysse et al. 
2015), as well as instrumental noise. Specifically, we consider 
contamination from the 88 GHz HGN(l-O) line, which we 
model in the same way as CO, with Ahcn = 1.7x 10“^® and 
bHCN = 5/3 (estimated based on Gao & Solomon (2004)). 
We assume the instrumental noise has a Gaussian PDF with 
zero mean and variance crjv. These contaminants can be 
added to the PDF by convolving P{T) for the original map 
with that of the contaminant. 

We base our model survey on the “Full” CO mapping 
experiment described in Table 2 of Li et al. (2015). This ex¬ 
periment surveys 6.25 square degrees between 2 = 2.4 and 
2.8, with CTjv = 5.8 p^. For simplicity, we assume that nei¬ 
ther the signal nor the foregrounds evolve significantly across 
the observed frequency range. We use information from all of 
the frequency channels to study the CO properties averaged 
over the full survey volume. 

Figure 1 shows the predicted one-point PDFs for the CO 
signal and the two contaminants. We use 50 logarithmically 
spaced bins between T = 2 and 1000 pK, neglecting intensi¬ 
ties outside of this range. The number of pixels within each 
bin is computed by integrating P{T) over the width of each 
bin. Below T ~ 1 pK, unphysical “ringing” effects come into 
the PDF due to the hard cutoff at Mmin- Since this regime 
would be noise dominated in a realistic experiment, the ef¬ 
fect on the CO constraints should be small. 


3 CO LUMINOSITY CONSTRAINTS 


In order to determine the constraining power of such an 
experiment, we perform a Fisher analysis on our calculated 
P{T) curves. We calculate the Fisher matrix 


1 dNkdNk 

'' dXi dXj ’ 

k = l 


(6) 


POO 

V{N)= / Pj^^{P)PvoissiN,P)dP, 
Jo 


where Nk is the predicted number of pixels in bin k and 
Xi = (Aco, 5co, Ahcn, bncN) are the four parameters 
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Figure 1. Predicted PDFs for intensity maps of CO, CO with 
HCN, and CO with instrumental noise {Thick curves: the full 



Figure 2. Predicted 95% confidence regions for CO mass- 
luminosity parameters Aqq and feco four cases: only CO 
with no noise or foregrounds (light blue), CO with foreground 
HCN (green), CO with instrumental noise (red) and CO with 
both noise and foreground (blue). 

we fit for. Note that we have assumed that the instrumen¬ 
tal noise is well understood a priori. The derivatives in the 
Fisher calculation are evaluated at the maximum fiducial 
parameter values given above. We have assumed Poisson 
uncertainties on the number of pixels in each bin, so the 
error on Nk is y/Nk- 

We can invert the Fisher matrix to obtain the covari¬ 
ance matrix, and marginalize over the foreground parame¬ 
ters to plot confidence regions for our CO parameters. Fig¬ 
ure 2 shows the 95% confidence ellipses for a map containing 
just CO, as well as for maps contaminated by HCN and in¬ 
strumental noise. With no foreground lines, the parameters 
Aco and feco are fairly degenerate, since to linear order in¬ 
creasing one or the other simply makes every halo brighter. 
Even with both foregrounds and noise included, the uncer¬ 
tainties are around the ~ 20% level in Aco and the ~ 1% 
level in hco, so relying on one-point statistics, these intensity 
maps can provide excellent constraints on the CO luminosity 
function. For comparison, Breysse et al. (2014) estimated an 
uncertainty on Aco of roughly an order of magnitude. 

Now that we have constraints on A and b, we can cal¬ 
culate constraints on other astrophysically interesting quan¬ 
tities, such as (Tco), which in addition to being a useful 
astrophysical quantity in its own right, is important for un¬ 
derstanding the power spectrum of an intensity map. The 
mean volume emissivity of CO emitters is simply 

POO J 

JCO = .400 (7) 

which can then be easily converted to brightness tempera¬ 
ture (Lidz et al. 2011). 


4 STAR FORMATION CONSTRAINTS 


The focus of this Letter is on high-redshift star formation, 
so we will now demonstrate how the P(D) analysis described 
above allows us to constrain SFRD from the CO luminos¬ 
ity function. Our mass-luminosity parameters A and b in 
Breysse et al. (2015) were originally derived from a set of 
empirical scaling relations discussed in Pullen et al. (2013) 
and Lidz et al. (2011). The CO and FIR luminosities of a 
galaxy have a well-known correlation which we write as 


^fir _ „ f L'co 
Lq ^VKkms-ipc^ 


-’fpiR 


( 8 ) 


where Cfir and Xfir are constants and Lco/Lq = 
4.6 X 10“®(I/co/K km s“^ pc^) (Wang et al. 2010; 
Carilli & Walter 2013). The FIR luminosity of a galaxy can 
then be related to its SFR through the Kenicutt relation 


SFR Lfir 

SeTF - 


(9) 


for some constant Csfr (Kennicutt 1998). The values for the 
above constants^ used in Breysse et al. (2015) are Cfir = 
1.35 X 10"®, Xfir = 5/3, and Csfr = 1.5 x 10"^°. 
Pullen et al. (2013) stated that the weakest part of this pro¬ 
cess is the relation between SFR and halo mass; here we 
assume they are related by a power law and use the above 
relations to write it in terms of A and b: 


SFR(M) = 9.8 X 10"^® f (1°) 


We can then integrate this over the mass function to obtain 
the mean SFRD ip in our survey, which turns out to be 0.12 
Mq/ yr/Mpc® for our fiducial model. 

Constraints on our original parameters Xi are straight¬ 
forward to convert to constraints on other parameters Yi. 
We simply need to multiply the Fisher matrix calculated 
in Equation (6) on both sides by the Jacobian matrix 
Jij = dXi/dYj. If Yi = {{Tco) I'P), we get the confidence 
regions plotted in Figure 3. The uncertainty on SFRD with 
noise and foregrounds included is on the order of ~ 10%. 
It should be noted that this is an optimistic calculation. 
When converting our Fisher matrix to SFRD and (Too), we 
assumed that the scaling relations used to connect Leo to 
SFR are well constrained. In practice, the uncertainties in 
these relations will reduce the constraining power of such a 
measurement. More galaxy observations will be required to 
reach the constraints shown here. However, it is clear from 
these calculations that the potential of a CO intensity map¬ 
ping experiment to constrain SFRD is quite high. 

Up until now we have been considering a survey cen¬ 
tered at a = 2.6. It is interesting to consider how the SFRD 
constraint varies with redshift, especially as the reliability 
of other methods diminishes with redshift. Figures 4 and 
5 illustrate this redshift dependence. For this calculation, 
we hold all of the instrumental parameters constant except 
the central observing frequency. Figure 4 shows the frac¬ 
tional uncertainty in SFRD as a function of redshift. Figure 
5 places the result with noise and foregrounds in context 
with results from the far ultraviolet (FUV) observations de¬ 
scribed in Madau & Dickinson (2014) and gamma-ray burst 
(GRB) observations described in Kistler et al. (2009). For 
Figure 5, we use a somewhat different form for p{z) than 


^ In the notation of Li et al. (2015), a = XpiR, /3 = log(CFiR), 
and (5 mf = 10^°Csfr 
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Figure 3. Predicted 95% confidence regions for mean CO bright¬ 
ness temperature (Tqq) and SFRD for the same four cases shown 
in Figure 2. 

we have used thus far. Our modeling, which is quite simplis¬ 
tic and only intended to serve as a proof-of-concept, gives 
a form for tplz) (dotted black line) which differs from the 
Madau & Dickinson (2014) fit (solid black line). To facili¬ 
tate comparison with the data (grey points), the blue curves 
in Figure 5 show the 1-cr uncertainty on SFRD calculated 
assuming the Madau & Dickinson (2014) form for ^/'(a). 

As mentioned above, there is theoretical uncertainty re¬ 
garding the relations in Equations (8) and (9). The magenta 
curves in Figures 4 and 5 show the effect of taking it into 
account. The dashed magenta curves in both figures assume 
a 10% uncertainty on Cfir and Csfr. Current results such 
as Carilli & Walter (2013) tend to give a constraint on Afir 
which is roughly an order of magnitude better than the con¬ 
straint on Cfir, similar to how in Figure 2 we obtain a 
much better fractional constraint on bco than Aco- These 
curves therefore assume a 1% uncertainty on Afir. These 
constraints are somewhat better than current values, but 
they are likely pessimistic compared to what will be avail¬ 
able once intensity mapping data are available. If there are 
resolved CO emitters in the same volume, one could cal¬ 
ibrate the CO-SFR relation and dramatically reduce this 
uncertainty. The solid magenta curve in Figure 4 shows a 
more optimistic scenario where the errors on Cfir and Csfr 
are 1% and the that on Afir is 0.1%. 

In order to test the model dependence of our results, 
we consider two aspects of more sophisticated CO emis¬ 
sion models which have been neglected thus far in our de¬ 
liberately simplified analysis. First, results such as those 
of Obreschkow et al. (2009) show that the CO luminosity 
function cuts off considerably earlier than the halo mass 
function. In Li et al. (2015), this effect comes into play be¬ 
cause the Lco{M) relation turns over at a certain mass. 
We take this into account by adding an exponential cut¬ 
off to our L{M) at a halo mass of M* = 2 x 10^^ solar 
masses, approximately where the Li et al. (2015) turnover 
appears. Secondly, semianalytic models such the one pre¬ 
sented in Lagos et al. (2012) and Popping et al. (2014) show 
that the CO luminosity of a galaxy depends on many param¬ 
eters besides its mass. Again following Li et al. (2015), we 
add lognormal scatter with a — 0.3 dex to our L{M) model 
to account for this, preserving some dependence on mass 
while allowing for fluctuations of other galaxy properties. 

The dotted curves in Figure 4 show the effects of these 
changes for the cases with and without instrumental noise. 
We can see that increasing the model complexity and sup¬ 
pressing the bright end of the PDF has relatively little effect 
at lower redshifts but noticeably worsens the constraints on 
SFR at high redshift, up to roughly a factor of 5 worse at 
z = 7. However, even at these high redshifts the overall 
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Figure 4. Predicted fractional I-ct uncertainties on SFRD for 
different target redshifts. The cyan, green, red, and blue curves 
show the same four cases as Figures 2 and 3. The dashed magenta 
curve shows the effect of adding a pessimistic 10% uncertainty on 
Cfir and Csfr and 1% uncertainty on Afir. The solid magenta 
curve shows the effect of reducing these model uncertainties to a 
more optimistic 1% and 0.1% respectively. The dotted cyan and 
red curves show the effect of adding an exponential cutoff and 
scatter to the power law L{M) model as discussed at the end of 
Section 4. These changes have little effect at low redshift, and 
reduce the constraint by a factor of ~ 5 at high redshift. 
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Figure 5. Comparison between our predicted SFRD constraints 
and existing FUV data from Madau &; Dickinson (2014) (grey 
points) and GRB data from Kistler et al. (2009) (orange points). 
The dash-dotted black line shows our simplistic fiducial ipiz) 
model, and the solid black curve shows the Madau & Dickinson 
(2014) fit to the data. The blue curves show the Tier uncertainty 
from CO intensity mapping calculated assuming the fitted ’4>{z), 
including foregrounds and noise, but excluding modeling uncer¬ 
tainty. The dashed magenta curves show the effect of adding 
the pessimistic 10% model uncertainty from the dashed magenta 
curve in Figure 4. 

uncertainty would still be dominated by the error on the 
CO-SFR conversion seen in the magenta curves. 

5 CONCLUSIONS 

From the above results we can get a good sense of the po¬ 
tential of this technique for constraining SFRD. In an ideal 
world with no foregrounds, noise, or modeling uncertainty, 
this 6.25 deg^ survey could constrain SFRD to ~ 1% at 
a ~ 3. This scenario is obviously unrealistic, but we could 
reduce the modeling uncertainty with more observations and 
reduce the noise with more sensitive instruments. It may 
be possible to reduce the foreground contamination as well. 
The masking procedure described in Breysse et al. (2015) 
will have little effect on the one-point statistics since it 
would simply involve ignoring all pixels greater than some 
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Tmax- However, as described for Lyo intensity mapping in 
Gong et al. (2014), it may be possible to reduce the fore¬ 
ground levels by using other data from the same area of the 
sky. 

Counterintuitively, the constraint on SFRD in Figure 4 
seems to improve with redshift when foregrounds and noise 
are neglected. This is because at higher redshifts, the same 
survey parameters correspond to a larger volume of space, 
which will include more sources. In our simple model, ^co 
and 6co do not evolve with redshift, so the only effect of 
moving to more distant sources comes from the evolution 
of the mass function. Eventually there will be few enough 
sources above Mmin that the constraint will worsen, but for 
our redshift range the volume increase outweighs this effect. 
When noise is included, the volume effect dominates at low 
redshift, but the amplitude of the noise relative to (Too) 
takes over at 2 ~ 5 or so. With foregrounds, the relative 
amplitude of CO and HCN decreases with redshift, starting 
to dominate over volume effects at around 2 ~ 6. The model 
uncertainties as we have considered them here are nearly 
redshift independent, so the constraint with these included 
depends weakly on redshift. 

As seen in the dotted curves in Figure 4, our constraints 
do not change drastically when we consider a more sophisti¬ 
cated model. The model used for the dotted curves contains 
significant scatter around the mean L{M) relation and adds 
a sharp cutoff to the luminosity of high mass halos, yet the 
change in fractional uncertainty on SFR is only order unity. 
Though even this more intricate model is still simpler than 
that of Li et al. (2015) or state of the art semianalytic mod¬ 
els, our qualitative conclusions should still hold. We leave a 
detailed quantitative prediction for future work. 

Currently there is substantial disagreement between 
high redshift SFRD’s measured using different methods, 
as illustrated by the roughly order-of-magnitude disagree¬ 
ment between the FUV and CRB data plotted in Figure 
5. (Kistler et al. 2009) found that this discrepancy can be 
explained if the FUV analysis underpredicts the quantity of 
low-luminosity galaxies. Since intensity mapping is far more 
sensitive to the fainter population than other methods it 
could provide a powerful means to resolve this discrepancy. 

An effect we have not taken into account here is contam¬ 
ination from continuum foregrounds. Though in principle 
these should be easy to clean out, it is unclear what form the 
residual contamination would take and how it would alter 
our results. Another effect worth considering is the depen¬ 
dence of these results on the beam size. It seems plausible 
that it would be easier to constrain the shape of the CO lu¬ 
minosity function with a smaller beam, since there would be 
fewer pixels with multiple galaxies. We leave consideration 
of both of these effects for future work. 

It may also be possible to constrain other quantities 
related to star formation with these intensity maps be¬ 
sides global SFRD. For example, because the convolutions 
in Equation (3) depend on the full range of possible lumi¬ 
nosities, the one-point PDF at all intensities depends on 
the chosen value of Mmin. Thus, even though the faintest 
pixels will be noise-dominated, it could be possible to deter¬ 
mine Mmin from our PDF. Another common unknown for 
high-redshift galaxies is what is the relative contribution to 
the overall intensity of high-mass vs. low-mass galaxies. This 
question is hard to answer with galaxy surveys, but it should 
be more tractable with intensity mapping. The methods de¬ 
scribed in this Letter should be readily applicable to these 
and many other problems of high-redshift star formation. 


From Figure 5, it is clear that CO intensity mapping 
has great potential for constraining the high-redshift SFRD. 
Though we are currently limited by modeling uncertainty, 
this should improve with better models and more observa¬ 
tions. The potential bound shown by the blue curves, though 
idealized, is remarkably strong despite using only a small 
survey area. In addition, the results shown above do not in¬ 
clude any additional constraints from the two-point statis¬ 
tics of a map. Though they neglect all possible line fore¬ 
grounds, Li et al. (2015) show that strong constraints on 
the CO luminosity function can be obtained from the power 
spectrum alone. Combining intensity mapping with exist¬ 
ing multiwavelength galaxy surveys should further improve 
the constraints. More work is necessary before CO intensity 
mapping can be used to accurately determine SFRs, but 
these results clearly demonstrate the value of this effort. 
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